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ABSTRACT 

Self-gravity becomes competitive as an angular momentum transport process in accretion 
discs at large radii, where the temperature is low enough that external irradiation likely con- 
tributes to the thermal balance. Irradiation is known to weaken the strength of disc self-gravity, 
and can suppress it entirely if the disc is maintained above the threshold for linear instability. 
However, its impact on the susceptibility of the disc to fragmentation is less clear. We use 
two-dimensional numerical simulations to investigate the evolution of self-gravitating discs 
as a function of the local cooling time and strength of irradiation. In the regime where the 
disc does not fragment, we show that local thermal equilibrium continues to determine the 
stress - which can be represented as an effective viscous a - out to very long cooling times, 
t c = 240£2~'. In this regime, it is also found that the power spectrum of the perturbations 
is uniquely set by this effective viscous a and not by the cooling rate. Fragmentation occurs 
for t c < j6 c ,i t Q _I , where jS cr i t is a weak function of the level of irradiation. We find that j6 c ,i t 
declines by approximately a factor of two, as irradiation is increased from zero up to the level 
where instability is almost quenched. The numerical results imply that irradiation cannot gen- 
erally avert fragmentation of self-gravitating discs at large radii; if other angular momentum 
transport sources are weak mass will build up until self-gravity sets in, and fragmentation will 
ensue. 
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1 INTRODUCTION 

Self-gravity may be important in the cool outer regions of both pro- 
tostellar and Active Galactic Nuclei (AGN) accretion discs, where 
the combined effects of pressure and shear cannot stabilize the flow 
against gravitational instab ility. The loca l linear stability of self- 
gravitating discs is simple ( lToomre|[l964h : it depends upon the pa- 
rameter, 



Q = 



nGY,' 



(1) 



where c, and £ are the sound speed and surface density of a low- 
mass disc (M disc <k M«), and k is the epicyclic frequency which, 
in a Keplerian disc, is the same as the angular velocity SI. The 
non-linear behaviour, however, is complex. A self-gravitating disc 
can either fragment into bound objects, or attain a "stable" self- 
gravitating state in which angular momentum transport results in 
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accretion and energy dissipation. Determining what sets the bound- 
ary between these out comes is import ant for understanding angular 



momentum transport ( Armitag el201lh and planet formation in pro- 
toplan etary disks i Boss 19971) . and star formation in discs around 



AGN l lBonnell & Ricell20oi) 

Since gravity is a long-range force angular momentum trans- 
port via self-gravity could in principle be a non-local pro- 
cess ( Balbus & Papaloizou 1999 ). Global numerical s imulations 
dLodato & Ricell2004 ICossins, Lodato & Clarke! l2009h . however, 
show that disc self-gravity can be approximated surprisingly well 
using a model in which turbulent heating associated with angular 
momentum transport is locally b alanced by radiat i ve losses. The 
use o f this local approximation jPaczvnski|[l978l : iLin & Pringld 
Il987l) greatly simplifies the description of the disc. In a local de- 
scription, thermal equilibrium implies that the stress, measured via 
the eq uivalent Shakura-Sunyaev a parameter dShakura & Sunvaevl 
1973), is inversely proportional to t he time scale o n which the 
disc can radiate its thermal energy dGammi iHoml). The stabil- 
ity of the disc against fragmentation also depends upon its cool- 
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ing time dShlosman & Begelmai] Il987l) . Writing i cool = fiQr 1 , 
two-dimensional local numerical simulations show that fragmen- 
tation of a self-luminous disc (one in which the only source of 
heat is furnishe d by accretion) ensues whenever p < /? crit = 3 
dGammid I200TI assuming a 2D adiabatic index y = 2). This 
minimum cooling time before fragmentation is equivalent to a 
maxi mum q crit ^0.1 that a stable s elf-gravitating disc can sus- 
tain ( Rice, Lodato & Armitagd 120051) . Consistent results for the 
location of the fragmentation boundary are obtained from two- 
dimensional global simulations, provided that care is taken to 
avo id prompt fragmentation during the ini tial growth of instabil- 
ity |Paardekooper, Ba ruteau & Merd 1201 lh . In three-dimensional 
global simulations the situation is less clear. There is evidence 
that such discs are somewhat less stabl e against fragment ation 
than their two-dimensional counterparts jMeru &Batdl201lh . al- 
though a critical cooling time does nonetheless appea r to exist 
dLodato & Clarkell201 it iRice. Forgan & Armitagell201 lh . We also 
note that strong temperature dependence of t he opacity can mod- 
ify th e fragmentation boundary substan tially Jjohnson & Gammiel 
120031 : ICossins. Lodato & Clarke] 120101) . This work deals exclu- 
sively with simpler temperature-independent disc cooling. 

Applying these results for self-luminous discs to proto- 
stellar and AGN discs, fragmentation i s predicted to occur for 
r Z 70 - 100 a u in protostellar discs Jivlatzner & Levinl 1 20051; 
Rafikovl 120051: rS tamatello s. Hubber & Whitworthl l2007t jBolevI 



20091 ; IClarkdl2009tlRafikovll2009tiRice. Mayo & Armitage|2010h . 
and at r > 0.06 (M BH /10 6 M n ) 1/3 pc for AGN jRafikovf 120091 ; 
Good man] 120031 ; ILevirj l2007h . At these large radii, however, the 
temperature of a self-luminous disc is small enough - about 10 K 
- that the thermal balance of a real disc in a star forming region 
or galactic nucleus may often be influenced by external irradiation 
(or, in the AGN case, by internal heating due to embedded stars). 
Here, we use numerical simulations to revisit the stability of self- 
gravitating discs, taking into account that irradiation as well as dis- 
sipation of accretion energy may contribute to the thermal balance. 
Our goal is to determine whether, in an irradiated disc, the fragmen- 
tation boundary is determined by a minimum cooling time or by a 
maximum a. These quantities, which are no longer equivalent in 
a disc subject to irradiation, are both physically important: [} mea- 
sures the time scale on which overdense clumps can radiate thermal 
energy and collapse further, while a is related to the amplitude and 
non-linearity of density perturbations. We also address the question 
of whether an annulus of the disc, initially assumed to be stabilized 
against self-gravity by irradiation, is able to viscously respond to 
an increase in surface density without enterin g the fragmentat ion 
regime. This work expands on earlier work bv lCai et"all J2008h in 
which irradiation from an envelope around a protoplanetary disc 
was shown to weaken the gravitational instability, and is com- 
pleme ntary to the analytic discussion in iKratter & Murrav-Clavl 
(|20 111) . 

The plan of this paper is as follows. In Section 2 we describe 
the local (shearing-sheet) simulations that we perform to investi- 
gate the evolution of irradiated self-gravitating discs, in Section 3 
we discuss the results of these simulations, and in Section 4 we 
discuss some implications of our results. 



2 METHOD 

To investigate the evolution of self-gravitating discs in the presence 
of external irradiation we use the Pencil Code, a finite difference 
code that uses sixth-order centered spatial derivatives and a third- 



order Runge-Kutt a time-stepping s cheme (see lBrandenburd fc003h 
for details). As in lGammid ( l200ll) . we use a "shearing sheet", or 
"local", model in which the disc dynamics is studied in a local 
Cartesian coordinate frame corotating with the angular velocity, fl, 
of the disc at some radius from the central star. In this coordinate 
frame, the unperturbed differential rotation of the disc manifests it- 
self as a parallel azimuthal flow with a constant velocity shear in 
the radial direction. We assume Keplerian rotation and hence use 
a shear parameter of q = 1.5 (such that the (/-component of the 
fluid velocity is u y = -qflx). The unperturbed background surface 
density, T,„, and two-dimensional pressure, P„, are assumed to be 
spatially constant. A Coriolis force is included to take into account 
the effects of the coordinate frame rotation. 

To use the "local" model we need to assume that t he disc is 
cool a nd hence thin (H/r = c s /(Cir) <k 1) and, as in iGammiel 
d200lh . we also assume that it is razor-thin. The continuity equa- 
tion and equations of motion in this model are 



— + V • (2k) - qflx— 
at dy 

du x du x 
— + (U ■ V)U X - qClx— 
at ay 







— + (u ■ V)Kj, - qClx— 



dt 



dy 



--— + 2£lu y -^-+f„ (2) 
l ax ox 

1 8P dd> 

where u(u x , u a ) is velocity relative to the background parallel shear 
flow « o (0, -qflx), P is the two-dimensional pressure, X is the sur- 
face density, <p is the gravitational potential of the gas sheet, and 
fvifvx, fnj) is a viscosity term that includes both a kinematic vis- 
cosity and a bulk viscosity for resolving shocks. Since this disc is 
razor-thin, Poisson's equation becomes 



: 4nGI.6(z). 



(3) 



This can be solved by Fourier transfoming (using the standard Fast 
Fourier Transform (FFT) technique) the surface density from the 
(x, y) plane to the (k x , k y ) plane giving 



<p(k x , k y , t) 



2nGL(k x ,k y ,t) 



(4) 



where (j>(k x ,k y ,t) and ~L(k x , ky,t) are the Fourier transforms of the 
gravitational potential and surface density, and k = ^jk^. + ky. It 
should be noted that in the shearing sheet approximation, the ra- 
dial wavenumber, k x , is no longer constant, but varies with time 
as k x (t) = k x (0) + qflkyt. For a two-dimensional domain of size 
L x L with resolution N x N, the only Fourier harmonics that will 
be present at time t will have wavenumbers k v = 2nn y IL and 
k x = 2mi x /L + q£lt'(2miy/L), where f = mod[t, l/(q£l\n y \)], and 
intege r number s n x an d n v lie in the range -N/2 < n x ,n y < N/2. 
As in Gammie (2001) we only use wavenumbers that satisfy k < 
nN/(Ly/2) which is the largest circular region in Fourier space 
that is always available, and ensures that the gravitational force is 
isotropic on small scales. 

The equation of state is essentially the standard 



P = (y- W, 



(5) 



where U is the two-dimensional internal energy per unit volume, 
and y is the two-dimensional adiabatic index. We use y = 1.6 
which can be mapped to a three-dimensional adiabatic index of be- 
tween ~ 1.6 and 1.9 dep ending on whet her the disc is strongly self- 
gravitating or not (e.g.. lGammi el d200lh ). For protostellar discs, a 
smaller value of y may be more appropriate. However, since the 
work here is not restricted only to protostellar discs, we have cho- 
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sen to use y = 1.6. Rather than solving for internal en ergy, the Pen- 
cil Co de actually solves for the specific entropy, s (see lBrandenburd 
(2003) for details). We use a cooling function similar to that used 
bv lGammid ( 1200 If) but modified by the inclusion of a heating term 
that represents external irradiation. When written in terms of spe- 
cific entropy, the combined heating and cooling function that is im- 
posed is 



ds 
dt 



y(y - 1)t c ' 



(6) 



where T is the temperature corresponding to a sound speed of c s 
and t c is constant cooling time, generally written as r c = /3fi~' 
with p a predefined constant. The sound speed, c so , is an effective 
minimum sound speed set by the background irradiation and which 
we generally express in terms of an effective Q ln where 

&n = -^T- (V) 

All the simulations presented here consider a rectangular box 



of size L x = L„ 



320 with a resolution of N x N where TV = 



1024, the same as the "standard" run in lGammiell200lh . The sound 
speed, or pressure, and surface density are initially constant and 
are set such that <2j„j t = 1. The initial perturbations are introduced 
through the velocity field which is perturbated from the background 
flow by a Gaussian noise component with a subsonic amplitude. 
Each simulation has a prescribed cooling time, /?, and an imposed 
level of external irradiation, gin - 




Figure 1. Quasi-steady surface density structure for/? =10 and Q m = 0. 



3 RESULTS 

3.1 Evolution without irradiation 

We consider initially simulations with no external irradiation 
(Gin- = 0) and vary the cooling time (as measured by ft) to es- 
tablish the fragmentation boundary and to investigate the quasi- 
steady nature of those simulations that don't fragment. A quasi- 
steady state is one in which Q settles to an approximately constant 
value, generally between 1 and 2, and remains at this value for 
many cooling times. In such a state, the instability is active but the 
system does not fragment and the instability acts to transport angu- 
lar momentum. The disc is regarded as having fragmented if very 
dense clumps form, with densities more than 2 orders of magni- 
tude greater than the average density, and survive for many cooling 
times. The presence of a bound clump also rapidly heats up the sur- 
rounding gas and so the Q value in a simulation that fragments also 
does not settle to an approximately constant value, but continues to 
rise. Fig. [T] shows the surface density structure from a simulation 
that has settled into a quasi-steady state. In this case /3 = 10 and 
<2in = 0. Fig. [2] shows the variation of Q against time illustrating 
that after an initial burst phase, Q settles into a quasi-steady state 
with Q ~ 1.8 that persists for many cooling times. Note that the 
saturated value of Q is almost a factor 2 larger than what is usually 
found in 3D simulations (e.g. Cossins et al 2009), where Q ~ 1.1. 
This is due to the stabilising effect of the finite disc thickness in 3D 



simulations dRomedl 19921 ; Mamatsashvili & Ricefe OlO). which di- 
lutes the effect of gravity by ~ ( 1 + kH), where k is the wavenumber 
of the perturbation. Since for the most unstable modes k ~ 1/H, we 
expect a reduction in the linear stability threshold for Q by approx- 
imately a factor 2, as observed. Fig. [3] shows the surface density 
structure for a simulation, with ft = 6 and Q m = 0, that has under- 
gone fragmentation. 

It can be shown that when a system settles into a quasi-steady 
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Figure 2. Q profile against time for fi = 10 and Q m = 0. 



state, the shear stress - wh ich we express in terms of an effective a 
- satisfies the relationship (Pringle 198ll; lGammiell200ll) 



(8) 



9y(7-l)T c n 



The shear stess can also be determined in each simulation using the 
Reynolds and gravitational stresses. The average Reynolds stress is 

(H xy ) = (Xu x u y ), (9) 

where u x and u v are, again, the velocity perturbations wit h respect 
to the background Keplerian flow. As described in detail in lGammiel 
(2001), the average gravitational shear stress can be determined in 
the Fourier domain using 



z 



\k\ 3 



(10) 



where the sum is over all Fourier components. The effective a is 
then 
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Figure 4. Comparison of the measured effective a values (crosses) deter- 
mined using equation jilt with that determined from the imposed cooling 
time (diamonds). 



Figure 3. Surface density structure for a simulation, with / 
0, that has undergone fragmentation. 



6 and Q m 



3(EcJ> 



(11) 



which can then be compared with the effective a determined using 
the imposed cooling time. 

We consider a series of simulations with /? v arying from 4 to 
240. We find that fragmentation occurs for/3 < 8.|G ammi e (2001) 
found fragmentation for /? < 3 and the reason for the difference 
is simply that we have considered a different specific heat ratio, 
y dRice. Lodato & Arm itage 2005$). For each simulation that does 
not fragment we calculate the effective a from the Reynolds and 
gravitational stresses using equation d I 1 b and compare it with that 
expected from the imposed cooling time (e.g., equation l[8j). This 
is illustrated in Fig. [4] in which the measured a values determined 
using the Reynolds and gravitational stresses are plotted as trian- 
gles, while that expected from the imposed cooling are plotted as 
diamonds. The averaged quantities that are needed to determine 
the effective a values are written out by the Pencil Code every 100 
iterations, corresponding to a time resolution of ~ Si At = 0.1, de- 
pending on the Courant condition. The measured a values are then 
determined by averaging over at least a few cooling times (gener- 
ally 5 or greater, except for the simulations with very long cooling 
times) starting once the system has settled into a quasi-steady state. 

There is a discrepancy between the measured and expected a 
values, with the mean of the measured values being higher than that 
expected. The lcr error bars, however, show that the values are at 
least consistent. The discrepancy is thought to be due to truncation 
errors. We have chosen to minimise the kin ematic viscosit y and as 
result some energy is lost at the grid scale JGammi 3 l200lh . Figure 
|4]does, however, show that the maximum stress that can be attained 
in a quasi-steady system is a ~ 0.06, consistent with earlier results 
dGammidl200ll ; lRice. Lodato & Armit age 2005). 



3.2 Evolution of irradiated discs 

In a quasi-steady state in the absence of external irradiation, the 
imposed cooling is balanced by heating due to the Reynolds and 



gravitational stresses. To consider the influence of external irradi- 
ation, we impose a combined heating and cooling function of the 
form shown in e quation <[6l). In a quasi-steady state the foll owing 
should then hold jGammij26bll ; lMamatsashvili & Ricell2009l) 

¥ + do £ f w *) = l Q + < G -» (12) 



Oct 



7(7 - IK- 

Using equation i ll It we can therefore show that, with external irra- 
diation included, 



9y(y- 1)Qt c 



<s>d 



(14) 



If we assume that (ScJ) ~ (£) (c,) 2 , we can rewrite equation d 1 4b 
as 

a 4 L Qlr 

~ <Wy-l)fiT c \ 



(15) 



where Q siit is the saturated Q value to which the quasi-steady sim- 
ulations settle. Since S and c, are probably correlated, the approx- 
imation made to go from Equation dl4t to Equation d 1 5 b may not 
be accurate, but at least gives us an approximate relationship be- 
tween a, Q m and Q sM . The form of equation d 1 5 b is also consistent 
with, a nd a generalisation of, the expression used by Lin & Pringle 
( 1990) to represent viscosity in a self-gravitating disc. 

To establish the influence of external irradiation, we consider 
various values of j} and for each value of f} we consider a number 
of different values of Q in . Fig. [5] compares the measured a values 
(determined using the Reynolds and gravitational stresses as de- 
scribed above) with the expected a values determined using equa- 
tion CO), plotted against Q m . for /3 = 7 (top left), fi = 8 (top right), 
and p = 9 (bottom). Again, we only include simulations that settle 
into a quasi-steady, self-gravitating state. The diamonds show the 
mean of the measured values and the error bars are 1 cr errors. The 
triangles show the expected value determined using Equation d 141 . 
Once again, the measured values are higher than those expected but 
do illustrate, very clearly, that as the level of irradiation increases, 
the instability weakens - as expected - and the value of a decreases. 
A quasi-steady a that exceeds the expected maximum of ~ 0.06 is 
also never measured. 

Fig. [6] summarizes our determination of the fragmentation 
boundary in the (<2i, r ,/6) plane. The plotted contours show lines of 
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Figure 5. Comparison of measured a values with that determined using the imposed cooling time, plotted against Q m for /? = 7 (top left), /} = 8 (top right), 
and/3 = 9 (bottom). 



constant a, determined using equation dl5l l with, using Figure [2] as 
a guide, (? s;lt = 1.9. Ideally we should be using Equation d 1 4b to 
set the contours, but we don't really know an appropriate value for 
(Sc^ or now it might vary with Q^, so Equation (JT3J provides a 
suitable approximation. The contour values, from bottom to top at 
Gte = 0, are a = 0.06, 0.05, 0.04, 0.03, 0.02 and 0.01. The symbols 
in Fig. [6] show all the simulations that have been performed, with 
the triangles being those that fragmented and the diamonds being 
those that settled into a quasi-steady state. The number next to each 
diamond is the measured a value for that simulation. 

Fig. |6] illustrates that for every cooling time there is a level 
of irradiation (<2hr) f° r which fragmentation does not occur. How- 
ever, the boundary between fragmentation or a quasi-steady, self- 
gravitating state does not appear to be strictly determined by a 
maximum value of a. If the boundary was determined largely by 
a maximum expected value of a, it should lie along one of the 
solid contours in Fig. [6] The actual boundary is illustrated by the 
dashed-line in Fig. [6] As [} decreases, the value of Q m required to 
halt fragmentation is larger than would be expected. The maximum 
quasi-steady a is also, consequently, smaller than expected. This 
suggests that the local cooling rate can influence fragmentation, al- 
though it appears that for every cooling time there is a Q ln that 
inhibits fragmentation. We would expect that for sufficiently large 
<2irr ( e -g-> vertical dash-dot line at Q ln = 1.95 in Fig. [6) there would 
be no instability growth as the disc would be maintained above the 
threshold for linear instability. Also, for Q m = 1.6, fragmentation 
occurs for both /? = 4 and /? = 5, while for Q m = 1.7 a quasi-steady 
state is achieved for both values of /?: this suggests that there may 
be a boundary region (i.e., between Q ln = 1.7 and Q m = 1.95) in 
which a quasi-steady, self-gravitating state is achieved for all values 
of p. However, it is also possible that for sufficiently small values of 
p there could be no value of Q m for which a quasi-steady (i.e., self- 
gravitating but non-fragmenting) state is achieved. If so, the sys- 
tem would either undergo fragmentation or be linearly stable. This 
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Figure 6. Figure showing contours of constant a determine using equation 
fBl with Q SM = l .9. The contour levels are a = 0.06, 0.05, 0.04, 0.03, 0.02, 
and 0.01. The symbols represent all the simulations, with the diamonds be- 
ing those that settled into a quasi-steady state and the triangles being those 
that fragmented. The number next to each diamond is the measured a value 
for that simulation. The dashed line is the apparent boundary between frag- 
mentation and a quasi-steady, self-gravitating state in the presence of irradi- 
ation. The dash-dot line illustrates the value of Q m above which we would 
expect systems to undergo no instability growth as they would be above the 
threshold for linear stability. 



could have implications for mass loading of systems with very short 
cooling times. If there are indeed cooling times for which a quasi- 
steady state cannot exist then it suggests that such systems would 
be unable to viscously respond to increases in surface density and 
would typically go from being linearly stable to fragmenting if the 
surface density increases sufficiently. The long growth timescales 
of the instabilty in systems with large values of Q- m is, however, 
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very long and so we don't attempt to determine here precisely what 
happens for very short cooling times. 

Ultimately, we have been able to robustly determine the frag- 
mentation properties of discs with < Q ln < 1,6, reaching values 
at the top end that are close to the limit where irradiation entirely 
quenches gravitational instability. In this range, we find that irra- 
diation modestly suppresses the critical value of /? below which 
fragmentation occurs, from a value of ^ 8 in the absence of 
irradiation down to a value y6 clit - 4 for the most strongly irradiated 
disc simulated. The maximum value of a that a self-gravitating disc 
can sustain similarly decreases as the strength of irradiation is in- 
creased 

ICossins. Lodato & Clarke! J2009h suggests that there is a rela- 
tionship between the perturbation amplitude (<5S/2) and the cooling 
parameter, p. In particular they find that the rms averaged perturba- 
tion amplitude, (<S£/£ avl , V obeys the following relationship 




If, as suggested above, the fragmentation boundary is determined 
primarily by the effective a rather than by the imposed cooling, 
P, we might expect that the perturbation amplitude would depend 
on a rather than p. Fig. [7] shows a plotted against {S'L/Y) for all 
the simulations we performed that did not fragment. The diamonds 
are for all the simulations with Q m = (i.e., Fig. [4]l while the 
other symbols are for the simulations with Q m t and corre- 
spond to the simulations illustrated in Fig. [5] although we have 
also added the p = 4, /? = 5 and p = 6 simulations that aren't 
included in Fig. [5] The values for ^<5E/X avg ^ were determined by 
averaging over at least 10 slices starting after the quasi-steady state 
was reached. The error bars are l<x errors. The curve in Fig. [7] is 
a oc ((t?£/£)) 2 and is essentially equiv alent to equation J 1 6b . taken 
from|c ossins. Lodato & Clarke! d2009h . Fig. [7] therefore confirms 
that the there is a relationship between the perturbation amplitude 
and the effective a, rather than the cooling time. That this holds 
for both irradiated and non-irradiated discs illustrates that it is uni- 
versal and is essentially independent of the condition of thermal 
equilibrium. 

The reason we have plotted a against <(5S/E>, rather than the 
other way around, is to highlight that the fundamental relationship 
is that the magnitude of the effective a depends only on the magni- 
tude of the density perturbations. If, as one might expect, fragmen- 
tation requires non-linear density perturbations ((c5£/£) > 1) this 
then implies that the maximum a that can be supplied by a quasi- 
steady system (i.e., one t hat does not fragment) is a ~ 0.0 6, con- 
sistent with earlier results dGammi el200ll ; iRice et al .l2003al) . How- 
ever, as illustrated above, it appears that for short cooling times 
fragmentation can occur even if the expected a is less than 0.06 
which suggests that large, localised perturbations in these systems 
can grow to form bound clumps even if the perturbation amplitude 
in the system as a whole is typically small. 

To further illustrate the universality between a and the pertur- 
bation amplitudes, we plot, in Fig. [8] the po wer spectrum of the 
perturbations for 4 of the simulations. As in iGammiel d200lh the 
x-axis is in units of kL/(2n) and since the power spectra all peak 
at kL/(2n) ~ 7, this illustrates that most of the power comes from 
wavelengths smaller than L/7 (i.e., the shear stress come primarily 
from wavelengths significantly smaller than the model size). Fig. 
[8] also shows that the power spectra all have the same basic form 
with, as expected, the amplitude increasing with increasing a. Fur- 
thermore, the two simulations with the approximately the same a 
values (but different values for p and gin) have almost identical 
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Figure 7. Figure showing the relationship between a and the rms averaged 
perturbation amplitude, (<SE/Z). The diamonds show the values for all the 
simulations with Qi rr = while the others symbols represents the simu- 
lations with Qj rr j= 0. The curve is a oc (((5X/X)) 2 illustrating that there 
is a relationship between the perturbation amplitude and a irrespective of 
whether irradiation is present or not. 
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Figure 8. Power spectrum for 4 of the simulations. The x-axis is in units 
of kL/(2n) so the peak at kL/ (2n) ~ 1 illustrates that most of the shear 
stress comes from wavelengths significantly smaller than the model size. 
The power spectra all have the same basic structure with, as expected, the 
amplitude increasing with increasing a. The two simulations with approxi- 
mately the same a values are almost identical illustrating that the structure 
of the turbulent state depends primarily on the effective value of a. 

power specta. This illustrates that the structure of the turbulent state 
does not depend on the cooling rate or the level of irradiation, but 
is determined only by the effective value of a. 



4 CONCLUSIONS 

For a self-gravitating accretion disc to fragment, gravitational in- 
stability must produce strong density perturbations that can cool 
quickly enough to collapse before shear and pressure forces dis- 
perse them. These requirements for fragmentation can be expressed 
via two dimensionless numbers: a maximum value of a, which 
fixes the amplitude of density fluctuations, and a minimum value 
of p, which measures the cooling time in terms of the local orbital 
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time scale. In the absence of external heating, local thermal equi- 
librium implies an exact equivalence between these descriptions, 
but this is no longer the case once irradiation becomes significant. 
In this paper, we have used local, two-dimensional simulations of 
self-gravitating discs to plot the fragmentation boundary in both 
isolated and irradiated discs. We find that the fragmentation bound- 
ary for irradiated discs does not lie at fixed values of either a or 
p. Irradiation necessarily weakens the strength of gravitational in- 
stability, but whether it stabilizes or destabilizes discs against frag- 
mentation depends upon the chosen metric. In terms of the cool- 
ing time scale, irradiation stabilizes discs: a strongly irradiated disc 
can avoid fragmentation for cooling times almost a factor of two 
smaller than the limit for a non-irradiated disc. In terms of the max- 
imum stress, however, irradiation has the opposite effect: the max- 
imum a sustainable in an irradiated disc is lower than that reached 
in the absence of irradiation. This is somewhat counterintuitive but 
is likely a consequence of the universal nature of the form of the 
power spectrum of the perturbations. Although irradiation weak- 
ens the instability, it is still possible for large, local perturbations 
to exist. If the cooling time is sufficiently short then these large 
perturbations could collapse to form fragments even if the typical 
perturbations are small (i.e., a c g < 0.06). 

For protoplanetary discs, the implication of our results is to 
reinforce the conclusion that the outer regions of protoplanetary 
discs are generally unstable to fragmentation at t he high accre- 
tion rates encountered soon after disc formation jClarkdl2009l ; 
ICossins. Lodato & Cla rke 2010). For a disc forming in a cold en- 
vironment (T ss 10 K), the rate of infall at radii r ~ 50 - 100 AU 
may well exceed the maximum accretion rate (M ~ 10~ 7 M Q yr ) 
that can be transported inward by the magnetorotational instabil- 
ity, even if the MRI i s active at these radii despite the effects of 
ambipolar diffusion teai & Stonel201 ]]). If so, continued infall will 
result in the build up of surface density until self-gravity sets in, 
at which point our results imply that irradiation will not save the 
disc from fragmentation. Irradiation will, of course, increase the 
surface density and temperature of the critically fragmenting disc, 
making it (even) more likely that the outcome of fra gmentation 
will be substellar objects rather than massive p lanets dRice et al.l 
2003bl: IStamatellos, Hubber & Whitworti|2 007). Discs forming in 
substantially warmer environments (T as 100 K) have a better 
chance of avoiding fragmentation, since the MRI in these systems 
may be able to transport all the infalling gas inward without locally 
exceeding the threshold surface density for gravitational instabil- 
ity. If all other aspects of star formation remained fixed, we would 
therefore expect to see more wide binaries in cool environments, 
and more young extended gas discs in warmer star forming climes. 

Finally, we caution that our results strictly apply only to the 
stability of isolated self-gravitating discs. Protoplanetary discs ap- 
proaching the threshold surface density for gravitational instability 
do so on account of infall, which need not be steady or axisymmet- 
ric. The dynamic al effects of i nfall appear to stabilize discs against 
fragm entation dKratter et al.l l2010l : lHarsono.~Al exander & Levinl 
l201lh . at least to some degree, but further work is needed to as- 
certain whether it is possible to avert fragmentation indefinitely at 
radii where isolated local disc models predict collapse. 
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